# Copyright (c) 2013, Chandler Armstrong (omni dot armstrong at gmail dot com)
# see LICENSE.txt for details

from collections import namedtuple
from itertools import combinations, islice
from random import shuffle

from networkx import DiGraph, Graph, topological_sort



MACHEPS = pow(2, -24)
E = MACHEPS * 10
N = 2**16 # max number of points in P


# utilities ###################################################################
def memoize(f):
    # memoization decorator for a function taking a single argument
    class memodict(dict):
        def __missing__(self, key):
            ret = self[key] = f(key)
            return ret 
    return memodict().__getitem__

def window(seq, n=2):
    "Returns a sliding window (of width n) over data from the iterable"
    "   s -> (s0,s1,...s[n-1]), (s1,s2,...,sn), ...                   "
    it = iter(seq)
    result = tuple(islice(it, n))
    if len(result) == n:
        yield result    
    for elem in it:
        result = result[1:] + (elem,)
        yield result


# geometry functions ##########################################################
@memoize
def det2((a, b, c, d)): return a*d - b*c

@memoize
def det3((a, b, c, d, e, f, g, h, i)):
    return a*e*i + b*f*g + c*d*h - c*e*g - b*d*i - a*f*h

@memoize
def circumcenter(((x_1, y_1), (x_2, y_2), (x_3, y_3))):
    x2_1, x2_2, x2_3 = x_1**2, x_2**2, x_3**2
    y2_1, y2_2, y2_3 = y_1**2, y_2**2, y_3**2
    a = float(det3(
        (x_1, y_1, 1,
         x_2, y_2, 1,
         x_3, y_3, 1)))
    b_x = -det3(
        (x2_1 + y2_1, y_1, 1,
         x2_2 + y2_2, y_2, 1,
         x2_3 + y2_3, y_3, 1))
    b_y = det3(
        (x2_1 + y2_1, x_1, 1,
         x2_2 + y2_2, x_2, 1,
         x2_3 + y2_3, x_3, 1))
    return -b_x / (2 * a), -b_y / (2 * a)

def midpoint((x1, y1), (x2, y2)): return (x1 + x2) / 2, (y1 + y2) / 2


# point location functions ####################################################
def sign(p, (a, b)):
    # return sign of p to a, b
    # a, b == origin, destination
    # left, right == +, -
    a_x, a_y, b_x, b_y = a.x - p.x, a.y - p.y, b.x - p.x, b.y - p.y
    sign = det2((a_x, a_y, b_x, b_y))
    if abs(sign) < E: return 0 # on edge
    if sign > 0: return 1      # left
    return -1                  # right

def pinE(p, E):
    # tests that p is left of all e in E
    # e in E is CCW  (positive concave)
    sides = 0
    for i, e in enumerate(E):
        sides += sign(p, e)
        if sides < i: return False
    return True

def plocate(q, p, D):
    q = list(q)
    while q:
        n = q.pop()
        if pinE(p, window(n+(n[0], ), 2)):
            successors = D.successors(n)
            if successors: q.extend(successors)
            else: return n


# edge legalization functions #################################################
def islegal(r, (i, j), k):
    # r is the new point
    # i, j is the edge to test
    # k is the vertex adjacent of r
    r_x, r_y, r2 = r.x - k.x, r.y - k.y, (r.x**2 - k.x**2) + (r.y**2 - k.y**2)
    i_x, i_y, i2 = i.x - k.x, i.y - k.y, (i.x**2 - k.x**2) + (i.y**2 - k.y**2)
    j_x, j_y, j2 = j.x - k.x, j.y - k.y, (j.x**2 - k.x**2) + (j.y**2 - k.y**2)
    return det3((r_x, r_y, r2, i_x, i_y, i2, j_x, j_y, j2)) < 0

def esearch(n_query, e, D):
    # n_query is the query triangle
    # e is an edge in n_query
    # returns the triange adjacent n_query at e
    q = D.predecessors(n_query)
    e = set(e)
    closed = set([n_query])
    while q:
        n = q.pop()
        if n in closed: continue
        successors = D.successors(n)
        if e.issubset(n) and (not successors): return n
        closed.add(n)
        q = D.predecessors(n) + q + successors # depth first search

def legalize_edge(r, e, D, sort=False):
    t0 = tuple((r,)+e)
    t1 = esearch(t0, e, D)
    if not t1: return # edge with no adjacent is always legal
    k = set(t1).difference(e).pop()
    i, j = e
    if sort:
        # OLD: sign does not handle p (arg 0) in J, so put k in e
        # WHY not just sign(k, (i, j)) == -1?
        if sign(j, (k, i)) == 1: i, j = j, i # CW of k would be CCW of r
    if not islegal(r, (i, j), k):
        # edge flip
        t2, t3 = tuple((r, i, k)), tuple((r, k, j)) # new triangles
        D.add_edges_from([(t0, t2), (t0, t3), (t1, t2), (t1, t3)])
        legalize_edge(r, (i, k), D)
        legalize_edge(r, (k, j), D)


# main ########################################################################
Point = namedtuple('Point', ['isreal', 'x', 'y'])

def DG(P, Q):
    # P and Q are lists of (x, y) tuples
    # P are points in D
    # Q are 'frame' points, an ordered quadrilateral that enframes all p in P
    Q = tuple(Q)
    D = DiGraph()
    P = [Point(True, x, y) for (x, y) in P]
    a, b, c, d = [Point(False, x, y) for (x, y) in Q]
    side = sign(c, (a, b))
    if side != sign(d, (a, b)): a, b, c = c, a, b # preserve order
    if side != 1: a, b = b, a # put a, b, in CCW order
    N0 = [tuple([a, b, c]), tuple([c, d, a])] # triangles are ordered triplets
    D.add_edges_from(zip([Q, Q], N0))
    shuffle(P)
    for r in P:
        t = plocate(N0, r, D)
        for e in combinations(t, 2):
            i, j = e
            side = sign(r, (i, j))
            if side == -1: e = e[::-1] # put i, j in CCW order
            D.add_edge(t, tuple((r,)+e))
            if side == 0: legalize_edge(r, e, D, sort=True) # edge case
            else: legalize_edge(r, e, D)
    return D

def VG(D):
    # build the voronoi graph from the delaunay graph
    V = DiGraph()
    for t in D.node:
        if D.successors(t): continue
        cc = Point(True, *circumcenter((p.x, p.y) for p in t))
        V.add_node(t)
        V_t = V.node[t]
        for e in window(t+(t[0],), 2): # sides in t
            s0, s1 = e
            s0_isreal, s1_isreal = s0.isreal, s1.isreal
            if not (s0_isreal or s1_isreal): continue # if not a
            a = esearch(t, e, D) # adjacent triangle
            if not pinE(cc, [e]): V.add_edge(a, t) # add all t's tests to a
            if not (s0_isreal and s1_isreal):
                b = (Point(True, 0, 0), Point(True, 0, 0))
            else:
                cc_a = Point(True, *circumcenter((p.x, p.y) for p in a))
                b = (cc, cc_a) # bisector of e
            if s0_isreal:
                if not s0 in V_t: V_t[s0] = set()
                V_t[s0].add(b[::-1]) # CCW order
            if s1_isreal:
                if not s1 in V_t: V_t[s1] = set()
                V_t[s1].add(b)
    return V

def vsearch(p, D, V):
    # locate point in V by finding point in D and transforming to V
    N0 = D[topological_sort(D)[0]]
    q = [plocate(N0, p, D)]
    n = None
    closed = set()
    for t in q:  # for triangle in q
        k = V.node[t].items()
        while k: # while nodes in triangle
            n_, tests = k.pop()
            if n_ in closed: continue
            if pinE(p, tests):
                n = n_
                q += V.successors(t)
                break
            else: closed.add(n_)
        closed.update(n_ for n_, tests in k) # close remaining nodes in t
    assert n # sanity check
    return n
